NACA TN 3661 


NATIONAL ADVISORY COMMITTEE 
FOR AERONAUTICS 


SELF SHIELDING IN RECTANGULAR AND CYLINDRICAL 

GEOMETRIES 

By Harold Schneider, Paul G. Saper , and Charles F. Kadow 

Lewis Flight Propulsion Laboratory 
Cleveland, Ohio 


TECHNICAL NOTE 3661 



Washington 

April 1956 


3897 


NATIONAL ADVISORY COMMITTEE FOR AERONAUTICS 


TECHNICAL NOTE 3661 

SELF SHIELDING IN RECTANGULAR AND CYLINDRICAL GEOMETRIES 
By Harold Schneider, Paul G. Saper, and Charles F. Kadow 


SUMMARY 

The steady-state diffusion P g and transport theory P 3 solutions 
of several self -shielding problems have been obtained for multiregion 
cells of rectangular and cylindrical geometries . These cells represent 
a given arrangement of moderator, fuel, and cladding) us ual ly considered 
to be homogeneous. The three -region cells studied were of slab geometry 
only and consisted of water, steel, and uranium regions. The two-region 
cells were of identical volumetric proportions and compositions, differ- 
ing only in the homogeneous dispersion of uranium in water. The dimen- 
sions of the equivalent two regions of the cylindrical cell were chosen 
so that the same homogeneous self-shielding factor of the corresponding 
slab cell was maintained. Additional results were calculated consider- 
ing the effects of the chemical binding in water on the flux distribu- 
tions and self -shielding factors. 

The neutrons were assumed to be monoenergetic) their distribution 
function assumed to be dependent on but one spatial coordinate, and the 
scattering in the center of mass system assumed to be spherically sym 
metrical. The P 2 approximation underestimated the magnitude of self- 

shielding effects relative to the P 3 approximation for all cases. Con 

sideration of chemical binding in water was unimportant in calculating 
self -shielding factors for the three-region cell but could not be neg- 
lected for the two-region cell where the water and fuel were homogene- 
ously mixed. An electromechanical differential analyzer was used to 
solve the P 2 and P 3 flux equations. 


INTRODUCTION 

In zero power criticality experiments, it is usual to mock up the 
reactor composition in the most expedient way. A fuel-element moderator 
assembly may be represented by fuel strips temporarily attached to 
plates of cladding material or by solutions of enriched uranium salt m 
liquid moderator. In analyzing these criticality experiments, it is 
important to know how a given fuel-element moderator assembly departs 
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from the homogeneous conditions. Inasmuch as the dimensions of repeti- 
tive cells in an assembly are generally very much less than the average 
mean free path, solutions of higher order than diffusion theory are re- 
quired to evaluate departures from homogeneity. 

The steady -state diffusion-theory solutions and the next higher 
order approximation of several self -shielding problems have been obtained 
for multiregion cells of rectangular and cylindrical geometry. The neu- 
trons have been assumed to be monoenergetic, and their distribution f unc- 
tion has been assumed to be dependent on but one spatial coordinate. 
Spherically symmetrical scattering in the center of mass system has also 
been assumed. The solutions of the diffusion theory Pg and the trans- 
port theory flux equations were obtained by a differential analyzer. 

(All symbols used herein are defined in appendix A.) 

The cases studied for rectangular geometry were (a) a three-region 
cell consisting of water, steel, and uranium, and (b) a two-region cell 
consisting of identical volumetric proportions and compositions except 
for the homogeneous dispersion of uranium in water. The case studied 
for cylindrical geometry was a two-region cell having the same homogene- 
ous self -shielding factor as that of case (b) . 

Additional results were obtained by approximating the effects of 
molecular binding in the water molecule by a method derived by A. 
Radkowsky in a classified publication. The thermal neutron flux distri- 
butions were used to compute the ratio R of total absorption in uran- 
ium to the total absorption in the cell. 

The detailed derivation of the flux equations in the Pj approxi- 
mation in rectangular geometry is given in reference 1. The equations 
for the P 3 approximation of cylindrical geometry are derived in detail 
in appendix B. These equations were also obtained by R. R. Haeffner in 
a somewhat different manner and presented in a classified publication. 


ANALYSIS 

Neutron Flux Equations and Boundary Conditions 
For Three -Region Slab Cell 

Assuming a constant isotropic source of thermal neutrons S A for 
water region A, the P,. approximation for this region leads to the 
equations (ref. l) 
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CD 

CO 

to 


Fj(x) + a 0 F Q (x) - S A = 0 
F^(x) + 2F^(x) + 2a 1 F 1 (x) = 0 
2F^(x) + 3F^(x) + 5a 2 F 2 (x) = 0 

3Fg (x ) + 7 a^F ^ (x ) - 0 


CD 


where F k is the k th (k = 0, 1, 2, 3) coefficient of the Legendre ex- 
pansion of the flux function. Physically, F Q is the total neutron flux 
and F-j^ is the net neutron current. The coefficients denoted by 

(k=0,l,2,3) are given by 

“K “ ** - 44, k 

where is the number of scattering nuclei per cubic centimeter of 

region A, 2 is the total macroscopic cross section, and o g ^ k is the 

j^th coefficient of the Legendre expansion of the scattering function. 

The values of Og k are calculated for spherically symmetric scattering 

in the center of mass system in appendix C. The value of is taken 

as 1 neutron per cubic centimeter per second. Identical sets of equa- 
tions hold in regions B and C with the term a replaced by b and 
Cj respect ively y ond. with sl zero source in these regions • 

Eliminating F x and F 3 from equations (l) and formally integrat- 
ing with respect to x 



with a similar set of equations in regions B and C. The equations in 
this form are suitable for solution by a differential analyzer. 
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The boundary conditions are 


Fq(0) 

> 

= 0 



F£(0) 

= 0 





► 

at the ex 

FqCt) 

= 0 


F^(r) 

I! 

O 

V 




F 0^A( a ^ 

= f o,b^ 


F 2 ^ A ( a ) 

= F 2,B^ a ^ 




F 6>> + P i, B W +2F 2,BW 


y at boundary a 


F 2,A^ _ F 2,B^ 


J 


(3a) 

(3b) 

(3c) 

(3d) 

(4a) 

(4b) 

(4c) 

(4d) 



(p) 


f o,b^ f o,c^ 

F 2,B^ = F 2,C^ 

+ ^.bCp) f 6,c^ + ^.qCp) 

b i c i 


at boundary p 


(5a) 

(5b) 

(5c) 


F 2,B^ a ) _ F 2,C^ 
b 3 c 3 

Equations (4) and (5) follow from continuity requirements on the 
F^ g at the interfaces a and p whereas equations (3) follow from the 

requirement that F 1 and F 3 vanish at 0 and r (ref. l) . 


(5d) 


By setting Fg = 0, equations (2) to (5) reduce to diffusion theory 

with the diffusion and absorption coefficients being given by D = —— 

3a-^ 

and &q, respectively. The specialization of equations (2) to (5) to 
two-region cases follows directly. 
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Specific slab configurations . - Figure 1 shows the rectangular 
cells studied and gives the specific dimensions used. The neutron cross 
sections are given in figure 2. The volume proportions of the repeated 
slab array used were steel: fuel: steel: water , 26.7:1:26.7:52.3 leading 
to symmetrical cell volume proportions of fuel: steel: water, 
0.50:26.7:26.15. The specific solutions obtained correspond to an en- 
riched fuel thickness of 3 mils (l.5 mils in the cell). 

The dimensions and nuclear constants listed in figures 1 and 2 and 
the values of a g k computed in appendix C were used to calculate the 

coefficients which were substituted into equations (2) to (5). For 
greater accuracy and convenience the strongly absorptive uranium region 
(2 = 33) was magnified 20 times by means of a linear transformation of 

SL 

the independent variable . 

Results for the chemically bound cases were obtained by the use of 
an effective hydrogen mass M and an average hydrogen transport cross 

section o\_ H , the latter having been derived by A. Radkowsky in order 
tr ^ 

to approximate the effects of chemical binding of the hydrogen. The 
quantity M is defined by 



where o c ^ is the average hydrogen scattering cross section averaged 

O 

over a Maxwellian distribution of neutron flux. 

TT TJ 

The values of M, and o g for water at room temperature are 

2.0, 31.4, and 46.7, respectively. For the unbound case, M was taken 
as 1, and as 25.6 barns from which it follows that a tr =8.55 

barns . 

From the flux distributions, the self -shielding factors R were 
calculated by 


R = 


total absorption by uranium 
total absorption in cell 


P ^ 

U uranium 


P fo(x)z“to + f Fj 

^ uranium ^ steel 


g teel MZ 




F U Z ,J V IJ + F steelj stee V teel 


' 0 a 


‘0 


steely 

a 
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where F 0 steel 

uranium region, 
negligible . 

The magnitude of self -shielding effects may be measured by devia- 
tions in R from the value of Rg = 0.712 calculated for the case of 

a constant flux throughout the cell, that is, F 0 U = F Q S , which is equiv- 
alent to assuming that all of the constituents of the cell are homo- 
geneously mixed. 

Analyzer procedure for slab solutions . - Equations (2) subject to 
the boundary conditions of equations (3) to (5) are solved by the dif- 
ferential analyzer in the following manner: 

The diffusion P 2 approximation results by setting Fg(x) = 0. 

From the boundary conditions, F^(o) must be zero. Since F Q (o) is un- 
known, it must be guessed in order to begin the solution at the boundary 
0. Using the initial F Q (o), the analyzer integrates the F Q equation 
until the interface a is reached. At the interface, F 0 is continu- 
ous. The discontinuous change of slope F^ B (a) that must be applied 

to initiate the solution for region B is calculated by equation 4(c) 
with F^ = 0. The F Q equation of region B is then integrated until 

the interface 0 is reached where the discontinuous change of slope 
*0 q ( P ' necessary to start the solution in region C is calculated by 

equation 5(c) with F^ = 0. The analyzer then proceeds to integrate 

the equations of region C until the outer boundary y is reached. If 
F q (t ) turns out 0, then the correct Fq(o) has been guessed. If 

F 6(r) i °> then another value F Q (o) is guessed until this criterion is 
satisfied. The final results are shown graphically in figure 2. 

After obtaining the diffusion answer, equations (2) were set up to 
be solved simultaneously. For the first trial, the diffusion theory 
F q ( 0) was used, and an initial value of F 2 (o) was guessed. The condi- 
tions of equations (4) and (5) were used to compute the discontinuous 
changes in the slope at the interfaces. Thereafter, both F Q (o) and 

F 2 (0) were guessed simultaneously until the conditions F^(r) = 0 and 

F 2 (r) = 0 were met. The total operating time for convergence turned 

out to be 8 to 16 hours for a two-region cell and 16 to 32 hours for a 
three -region cell. 


is the average flux in steel, V U is the volume of the 
The absorption in water is neglected since Z H 2° is 

cl 
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Equivalent Two-Region Cylindrical Cell 

As shown in figure 3, the radii of the two -region cylindrical cell 
were chosen such that the volume of water plus uranium and the volume of 
steel were equal to those in a square cylindrical cell formed from the 
dimensions used in the previous two-region slab cell. With an outer 
square boundary, the neutron flux will be spatially dependent on an azi- 
muthal angle in the cross section of the cylinder in addition to radial 
dependence. This angular dependence is not present when the outer square 
boundary is approximated by a circular one. 

From figure 3, the volume of steel per unit depth for the rectangu- 
lar cell is 0.32773 cubic centimeters, and the volume of water plus ur- 
anium is 0.32098 cubic centimeters for a total cell. volume of 0.64871 
cubic centimeters . This means that for the cylindrical cell 

na 2 = 0.32773 


or 


a = 0.3230 cm 

which is the radius of the steel region, and 

rtb 2 = 0.64871 


b = 0.4544 cm 


the radius of the cell. The self -shielding factor R H then remains the 
same as in the slab cell (0.712) for the case of a constant flux through- 
out the cell. 


The cylinders are assumed to be- infinitely long and the transport 
flux depends spatially on the radial distance r from the z -axis, and 
on two angular coordinates a and 5 that determine the direction o 
the neutron velocity vector (see fig. 4). The azimuthal ® ° 

the velocity vector is measured from any fixed radial line. in the cros. 
section of the cylinder. The polar angle £ of the velocity vector is 
measured from the z-axis. The transport flux can be written as 
F(r,p,a) where p = cos £. 

In appendix B, the transport flux is expanded in spherical harmon- 
ics to give 

CO oo 

w-r cl ') - Y y (n-m) j. 2n+l ^ ( r ) CO s m a P"j(p) ( BQ ) 

F(r,p,a ) - L_j ( n+m ) ; 2 m,n' 7 n 

n=0 m=0 
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From the orthogonality of spherical harmonics, 

/ TC 1 

( 2lt8 j,0 + rt6 j^ F j,k( r ) = f F(r,p^a)cos ja P*^(n)dp da 

where the 6's are defined as 



When j=k=0, the right hand side is the transport flux integrated over 
all directions. This is the total flux $ 0 Q (r) . Hence 

= ^ rtF o,0^ r ^‘ ^kh J=k=l, the right side is the projection of 

the vector flux F(r,p,a) on the radial axis integrated over all direc- 



Assuming a constant isotropic source of neutrons in each volume 
element of the water-uranium region, in the P3 approximation (see 

appendix B), the following set of equations is obtained: 


F ' n + 2b F = 2S 

1,1 r 0 0,0 0,0 


3 

F ' _ + 

1,3 r 


1,1 - ^ + 10 Vo ,2 " 0 


?J 6,0- 2F 0,2 + 2 F 2,2 + ~r~ + 3b l F l,l " 0 


(B9a) 

(B9b) 

(B9c) 


1 r 2 2 

12F ' - - F' + 7b F = 0 

0,2 2 2,2 r + 3 1,3 


(B9d) 


F 3,3 + 3 -¥ + 12 


( F i,i - V) - 2 (5, 


F 

3 r 




15 ( F 2,2 - 2 - - 0 


3 3,3 


(B9e) 
(B9f ) 
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where b is defined in terms of nuclear constants by 


co 

6 

o 


b k - No - N g a s ^ k 


k = 0,1,2, 3 


Transforming to fluxes k and the neutron source s Q ^ 0 by 


$ , = 2*F^ , 

0,k 0,k 

# . , = rtF . , 
j,k 


S 0,0 “ 2 ltS o ,0 


k = 0,2 

j = 1,2,3 
k = 1,2,3 


and after eliminating F lj(1 , F^ 3 , and F 3 ^ 3 from equations (B9) and 
formally integrating with respect to r yields 


- 3^[*o,o + f^\ + b o f *o,o = So ,o r + 357 [' *6,2 + + if 

- ^[*0,2 + /^] + 0,2 “ S 0,0 r - b o/*0,0 - Tibi [*2,2 + 

- (5 + %][*2,2 + Kfl] + 10b 2/*2,2 

-ME 


dr 


V(6) 


The first of equations (6) will hereafter be called the $q, 0 equation. 
The second will be called the i>Q g equation and the third the ^ 2,2 
equation. 

An identical set of equations holds in region A (see fig. 3), ex- 
cept that Sq q = 0 and the terms containing b are replaced by terms 

containing a . 

Boundary conditions. - Referring again to figure 3 the following 
boundary conditions must hold: 


F(b,p,a) = F(b,-p,a+rt) 

at r = b 

(7a) 

F A (a,p,a) = F B (a,p,a) 

at r = a 

(7b) 
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F(0, n,^) = F^^ia^ctg) at r = 0 (7c) 

and arbitrary 

F(0,p,co) finite 

Equation (7a) is a statement of the fact that no net flow of neu- 
trons occurs across the outer boundary of the cell in any arbitrary 
direction. Though this condition is strictly true for cylindrical cells 
with square outer boundaries, it can only be approximately true for 
cells with circular outer boundaries. 


Equation (7b) states that the transport flux in an arbitrary direc- 
tion at the interface must be continuous across that interface. 


Equation (7c) is necessary because at r = 0 for a given polar 
angle the flux must be independent of a. This leads to the same con- 
ditions on the F m ^ n (r) as the requirement F(0,p,a) be finite. 

Equations (7a) to (7c) lead to the conditions 


- F 3,3( b > ■ 0 

F® (a) = F® (a) 
m>n'' m,n' 

- F l, 3 <°) - F 3 , 3 C°) ‘ F 2 , 2 <°> - °J 

which imply 


►m,n«0,l,2,3 


<i> f = -5$ 1 

* 0,0 ° * 0, 2 

$ 

0,2 3 b 


at r = b (outer boundary) 


2$ 


$ 2,2 b 


2^2 


(8a) 

(8b) 

(8c) 


a l ( a l 

a> ' . = — * ' + _ _ 

*0,0,A b x *0,0, B \b x 


bJV 2 0 


, \ i f 2a^ \ 

3,2,bJ + 3S \^~ + bj" ' 5 J*2,2,B 


a 3 1 

d) * = — + — 

*0,2, A b 3 0,2,B 3a 


( x ' 


B 


(9a) 

(9b) 
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a 3 2 A f3\ 

,2,A = bJ*2,2,B + a ^ ' b 3 J f 2,2,B 

at r = a (interface) 


a A - 6 B 

* 0,0 " 0,0 

* 0,2 = ^ 0,2 

* 2,2 ' * 1,2 


(9c) 

(9d) 

(9e) 

(9f) 


<|> ' ^ == 0 

0,0 

4 > * =0 

0,2 

$ ’ =0 

2,2 

4 s* 0 

* 2,2 


at r = 0 (origin) 


(lOa) 

(10b) 

(lOc) 

(lOd) 


ruffpr-ftntial analyzer procedure for solv i ng cylindrical fl u^c 


equations . - A neutron source of s Q ^ 0 
(second) is chosen in region B with a zero source in region A. ^s 
$ ^ 2 are set equal to zero in equations (6) to (9), and the 

resulting cylindrical diffusion equation is solved for *0,0’ 

The diffusion equation for region A and the one for region B are 
solved a m^er similar to that described in the section Analyzer 
Procedure for Slab Solutions. This procedure gives the graphical 2 

diffusion theory answer. 


= 1 


neutron/ (cubic centimet er ) 


Having the Pr 


solution, the procedure is to successively solve 
. * * * $ , $ A , and so forth, equations by using 

the *0,2' * 2 , 2 ' %0' *0,2' 2,2' 0 0' ^ ^ equation, 

the boundary conditions > equations v ) \ )> 0^2 

the third conditions for the * 2#2 equation, and the first condition 
for the q equation. 

Unlike the rectangular case, inhomogeneous boundary abactions occur 

at r =17 but the solution is still initiated by gues, sing < Kb) . The 

discontinuous change in slope at the taterface is c^culated^^ the 
responding $ equation of equation (8). As in tne recuangux ' 

homogeneous conditions on the slopes of the terms ♦ occur at 0. 

However, the additional requirement * 2 , 2 C°) = 0 must a 
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For the cylindrical case, it was not possible to set the three 
equations of equation (6) into the analyzer at the same time. There- 
fore, the following procedure was used: 


(1) ^2,2 seb = 0 • The 0 dr obtained from the diffusion 

§ 0 o equation is fed into the <f> „ equation at the proper time. A 

/ . . (l) 

first approximation ^ denoted by <f>Q g is obtained when the value 
assumed for $ 0 2 (b) leads to ^ 2 (o) » 0. 

(2) The functions (#' _ - P — ^ j obtained from and 

(v -PY 


obtained from the diffusion <t> 

t-g g cquo.oj.un. o ^*2^2 
t0 *2,2^ = °> 3111 $ 2 2^ = ° 


„ equation. A is obtained when the 2 (b) assumed leads 


0,0 
2,2 


were fed into the 


(3) The functions 


f m* 2 A ^ ^ ( i 

[$2 2 + 3 J — ^—1 obtained from $); 




r 

(1) 

2 


obtained from and 

0 , 2 ’ 

were fed into the <f> Q Q equation. 


Equation is obtained when the 0 (b) assumed leads to 


% 0 (°> " °' 


/* 


( 4 ) The functions (<f>^ ^ + 3 


.tl) 




obtained from $ 
were fed into the <l> 


( 1 ) 


and 


0,2 


2,2 

equation. The 


n n dr obtained from . 

Vjv . . 0,0 

( 2 ) 

term $ 0^2 i s obtained when the g(b) assumed leads to 4*^ g(o) = 0. 
Steps 2, 3, and 4 were repeated until the $ 0 Q equation solution pro- 
duced a neutron conservation of approximately one percent. In solving 
( 2 ) & 
for ^ 2 , 2 ^ successive values of $ 2 . 2’ interpolation of the function 


k. -f¥) 


appearing in the 3v> ^ equation was necessary. The 
curve to be used for this function had the same shape but it was inter- 


mediate between 


k. -F¥ 

k (2) / 


.(1) 


and 


k. -f^V - k- 


(*6 , 2 ~S * 8X1(1 between 

\( 3 ) 

and so forth. The amount 
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of Interpolation of the function ($q^ 2 


-y» 


decreased as the error 


in conservation of the 4 >q q equation decreased, until the last curve 

of this function used was practically that obtained from the preceeding 
a , equation. For the two cylindrical problems solved, this procedure 

required about 40 hours each. 


RESULTS AND DISCUSSION 
Rectangular Case 

The flux distributions for the slab cells studied are shown in fig- 
ure 2, and table I presents the numerical values of the self -shielding 
factors. The Pg approximation differs from the P 3 approximation in 

calculating self -shielding effects for the examples presented. 

The magnitude of self -shielding effects are measured by deviations 
in R from the value for R H of 0.712 calculated by assuming a constant 

ppux throughout the cell. For the two-region cell, the deviation from 
R h was 0.0107 and 0.0150, respectively, for P 3 , unbound and bound, and 

0.0048 and 0.0080 for Pg, unbound and bound. For the three-region cell, 
the deviations were -0.0225 and -0.0224 for P 3 , unbound and bound, and 
-0.0121 for Pg, unbound and bound. 

The deviations are positive for the two-region cell because the 
presence of the source term in the water-uranium region leads to a higher 
flux and, consequently, more absorption in uranium relative to the case 
where the flux was constant throughout the slab cell. In the three - 
region cell, the neutron flux is considerably attenuated by transition 
through the water and steel regions before absorptions occur in uranium. 
Hence, a negative deviation from Rg is obtained. 

A consideration of chemical binding does not affect values of R 
for the three -region cell despite an increased flux in water resulting 
from larger average hydrogen scattering cross sections for the bound 
cases The flux is higher in the water region for the bound cases be- 
cause 'neutrons undergo more scattering collisions per second which tends 
to confine them there and build up the flux. This effect is enhanced to 
some extent by the larger mass number M for the bound case which tends 
to produce more nearly isotropic scattering that attenuates the stream- 
ing' of neutrons toward the steel region. Thus a higher R is to be ex- 
pected for the bound cases. In the two-region cell, chemical bin ing 
leads to a larger value of R by virtue of the increased flux in the 
water -uranium region relative to the unbound case. 
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Because the total neutron flux is symmetrical about the extreme 
boundaries of the cell, no leakage occurs from the cell. Thus, for 
steady state, the total number of neutrons absorbed per second in the 
cell must equal the total number produced per second. This criterion 
was applied to each of the flux curves of figure 2 as a check on the 
accuracy of the analyzer solutions. As can be seen from table I, con- 
servation of neutrons in the slab case is satisfied within about 2 per- 
cent for the three-region case and a l/2 percent for the two-region 
case; the value of R was found to be insensitive to such small devia- 
tions in conservation. 


Cylindrical Case 

Unbound . - The P 3 approximation <f> 0 Q curve and the diffusion 
$0,o curve shown in figure 5 resulted from solving equations (6) to (9) 

for the equivalent unbound two-region cylindrical cell problem. The 
cylindrical P 3 $ Q Q com P ares with the P 3 approximation flux F Q (x) 

for the slab cell in the following manner: 

The values of $ Q _,o at tile outer boundary, interface, and origin 

are 1.184, 1.111, and 1.023, respectively, for the cylindrical, cell. 

The corresponding values of F Q are 1.071, 1.105, and 1.160. The flux 

in the steel region of the cylinder is lower implying a higher flux in 
the water -uranium region to satisfy conservation. 

Table II gives the results for the cylindrical cell. Unlike table 
I, the average production and absorption are taken per cubic centimeter- 
second of material. The average absorption per cubic centimeter -second 
of region A, for example, is calculated by 



After calculating the average absorption per cubic centimeter-second in 
region B, it follows that R » 0.7198 for the P 2 approximation and 

R « 0.7336 for the P 3 approximation. Thus, more absorption occurs 
in the water -uranium region of the cylindrical cell than in the same 
region of the slab cell. 
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As in the slab-cell case, conservation of neutrons must be satis- 
fied. The diffusion <f> 0 ^ 0 and the P 3 <^ 0 were approximately 2.0 and 

1.2 percent, respectively, too small conservation wise. As for the rec- 
tangular case, values of R were found to be insensitive to small devia- 
tions in conservation. 

Bound case. - The P 2 and P 3 curves for the bound case are shown 
in figure 6. As was explained for slab cells, a higher flux in the water - 
uranium region, leading to a larger value of R relative to the unboun 
case is to be expected when the effects of molecular bidding in water are 
considered. This is verified from table II, where R = 0.7213 for 2 
and R = 0.7355 for P 3 - Conservation is 2 percent and 0.7 percent lack- 
ing for these curves . 

Values of $q q at the outer boundary, interface, and origin are 

1.200, 1.111, and 1.018, respectively. The corresponding values for 
F (x) in the P 3 approximation for the rectangular cell are . , 

1.092, and 1.180, respectively. 


CONCLUDING REMARKS 

This study shows that the P 3 approximation of transport theory as 

compared with the diffusion theory approximation gives significant dif- 
ferences for the flux distributions and self -shielding factors of xml 
region cells of rectangular and cylindrical geometry. rela _ 

approximation underestimated the magnitude of self -shielding effects rel 
tive to the P 3 approximation for all cases. 

Also the chemical binding of the hydrogen in the water molecule was 
found to be important for the two-region cells which considered the ur- 
S to E homogeneously mixed with the water, 

. . Q . opi f- shield ine factor were affected. For tne -cnree regiuu 

cell where the uranium was considered to be a region, the .f^te ^ tlT " 
tog effect did not change the self -shielding factor but did alter the 

flux distributions somewhat. 

An electromechanical differential analyzer was found to be suitable 
for solving the P 2 and P 3 flux equations. 


Lewis Flight Propulsion Laboratory 

National Advisory Committee for Aeronautics 
Cleveland, Ohio, January 17, 1956 
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APPENDIX A 
SYMBOLS 

A,B,C 

used as either subscripts or superscripts to denote the var- 
ious regions of a cell 

a,b 

radii of the cylindrical cell 

vv c t 

defined as the difference (z - Ng ^ ) 

D 

the neutron diffusion coefficient defined by D = -i— 

J Sa-L 


the coefficients of the expansion of the flux function in 
cylindrical geometry 


the k th coefficient of the Legendre expansion of the flux 

function in rectangular geometry " 

F(r,n,a) 

definition of the transport flux in cylindrical geometry 

H 

chemical symbol for hydrogen 

* <J ^ 

i^j^k 

unit vectors along x,y,z-axes 

k 

0,1, 2,3 

M 

atomic mass 

% 

number of scattering nuclei per cubic centimeter of a region 

pm 

n 

the associated Legendre polynomial of indices n,m 

R 

self -shielding factor defined as the ratio of absorptions in 
uranium to the total number of absorptions in a cell 

r h 

homogeneous self -shielding factor, R^ = 0.7120 

r 

radial distance from the z-axis 

r 

unit radial vector 

S 

constant isotropic source of thermal neutrons, -p- 

7 2n 

U 

chemical symbol for uranium 
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V 

x 

0,a,|3,& 


a 

8 

0 

0 

c 

M- 

v 


K 

l 



$ 

«p 

£ 


volume 

variable distance used for the rectangular cell 

values of x at the various boundaries of the rectangular 
cell 

azimuthal angle of the velocity vector measured from the 
z-axis 

symbol for the Kronecker delta 

scattering angle in the laboratory frame of reference 
scattering angle in the center of mass system 
defined by |i = cos K 
defined by v = cos 0 C 

polar angle of velocity vector measured from the z-axis 
the total macroscopic neutron cross section 
macroscopic neutron absorption cross section 
macroscopic neutron scattering cross section 
microscopic neutron scattering cross section 

the k th coefficient of the Legendre expansion of the scatter- 
ing function 

microscopic neutron transport cross section in the laboratory 
fr am e of reference 

transport flux 

angle between radius vector r and x-axis 

azimuthal angle made by a neutron direction after collision 
with the neutron direction before collision 

unit vector in the direction of neutron velocity 
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APPENDIX B 

THE P 3 APPROXIMATION FOR CYLINDRICAL GEOMETRY 

The general steady-state monoenergetic Boltzman transport equation 
stating the conservation of neutron flux F for an element of volume is 
(ref. 2) 

2*VF(?,2) + Na(2)F(r,2) = P N s a s (2, 2* )F(r,2' )d2' + s(r,2) (Bl) 

vq' 

t 

d2» 

N 


a 

probability per centimeter that a neutron traveling in the 
direction 2 ' undergoes at r a scattering collision 
into a unit solid angle about direction 2 

The first term represents the net number of neutrons with direction 
2 at r leaking out through the faces per cubic centimeter per second. 

The second term is the loss of neutrons of direction 2 at r per cubic 
centimeter per second due to absorption or scattering collisions. The 
third term is the number of neutrons per cubic centimeter per second 
traveling in the direction that undergoes a scattering collision 

into the direction 2 at r . The last term is the number of neutrons 
produced per cubic centimeter per second with direction 2 at r. 

The following simplifying assumptions are made: 

(1) The neutron flux has cylindrical symmetry, that is, F is in- 
dependent of and z, the spatial azimuthal and axial coordinates for 

directions of 2 are fixed relative to the direction of r. 

(2) The medium is isotropic, which implies a (2) = o, a constant, 
and the scattering is a function only of the angle between the initial 
direction of the neutron velocity and the final direction. 


unit vector in the direction of neutron velocity 
element of solid angle about direction 2 1 
number of nuclei per cubic centimeter 
number of scattering nuclei per cubic centimeter 
total microscopic neutron cross section 
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(3) The scattering is spherically symmetrical in the center of mass 
system. 

Referring to figure 4 in which the direction of the neutron veloc- 
ity vector at r is given by polar angle K and an azimuthal angle a 

in the plane of the cross section, 

2 = "i sin 5 cos (a + <p) + "j sin £ sin(a + <P) + k cos K 

= r sin K cos a + <P sin € sin a + k cos K (B2) 


where r is a unit radial vector and (p is a unit vector perpendicular 
to r in the plane of the cross section. Then 

-*• . -*■ 3iF <P dF , &f\ 

2 -VF = (r sin K cos a + <P sin 5 sin a + k cos ' l r '^r + rc*p + K '^zJ 

ri 1 ^ 

=s sin £ cos a F(r,£,a) + — sin £> sin a F(r^^-^a) 4* 

cos € ^ F(r,£,a) 


In computing sin F(r,£,a), 2 must be held fixed as <P 

changes. (With 2 held fixed, K and- a do not vary with r, or z, 
as either variable alone is changed.) This requires that E, and a 
be functions of <P despite no explicit dependence of F on <p. 
Therefore 

2-7F = sin £ cos a || + i sin £ sin a||^F(r,£,a) || + ^ F(r,£,a) ^ 

(B3) 

since cos £ ■J” F = 0. 

The condition that 2 remain fixed as <p changes requires that £ 
and <p + a, the angles made by 2 and the fixed z- and x-direction, 
respectively, be constant so that 


(B4) 
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Letting p - cos E } substituting (B4) into (B3) and noting that 
assumption (2) implies 

a s (a,S2') = — a s (p 0 ) 

— ► . 

Pq = cos^SZ'J = cos E cos E' + sin E sin E’ cos(a - a') 
equation (Bl) becomes 

1/2 


/■-, 2\ ' f ciF sin a dF~\ „/■ » 

(i - P ) <cos a ^ — V + N a F(r,p,a) = 

n 1 n * % 

J 1 J 2ir a s^0^ F *‘ r ^ , -' a, ) da ' d ^' + S(r,p,a) 

a s (p 0 ) is expanded in a series of Legendre polynomials to give 


(B5) 


°s^o) " Z 


1=0 


2Z+1 


2 S,Z Z VH 0 


°s 7 F 7^n) 


where 


a s,z ~ J~' 1 a S^0^ P Z^0^ d ^c 


It is to be noted that Cg q is the total scattering cross section and 
o s i is the total scattering cross section times the average cosine of 
the scattering angle. 

By the addition theorem for Legendre polynomials, 


Vs , , _ ± 

2rt Uty “ 2it 


s 


E 

z=o 




oo l 


Z Z {z+m}: (2Ul)a S,Z P Z^^') cos M(a-a’) 


z=o ]yt=i 
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which gives the variation of the scattering function in terms of the 
desired variables. Equation (B5) becomes 


(1 - H 2 ) 


{ dF sin a dF 
cos a ^ — 


+ No F(r,p,a) = s(r,|i,a) + 


X) [z+m H C 2i+i )° s z p i^) p i^') cos m(<x " a, Y da ' d * A ' ^ 

J=OM * 1 


Symmetry demands that 

F(r,n,a) - F(r,|a,-a) 

F(r,p,a) =* F(r, -p,a) 

The transport flux is expanded in a set of spherical harmonics to 


(B7) 


give 


oo ao 

P(r,„a) - E C ° S “ Vn 5 0 < B8 > 


n=sO m=0 


m > n 


Harmonic terms containing sin a do not appear in (® 8 ) because oP t ^ e 
first of conditions (B7). The second of conditions (B7 ) implies that 

-p fr) = 0 for either subscript even and the other odd since 
m,n v ' 

P m (-M’) = (-l) n+m F^Cm-) • In the P 3 approximation, F m ^ n C r ) is neglected 

when n> 3. The function S(r,p,a) is likewise ejqpanded in spherical 
harmonics (for isotropic sources only the m = 0 terms are present;. 

Multiplying both sides of equation (B8) by cos ja P^. (p)dp da 
and integrating over all p and. a gives 

(* v + * = ££ Ftr ’“' a)cos Ja P * ^ da 
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j , k - 1,2,3 in the approximation. The ordinary Kronecker delta 

5^3 is defined as 



m =4 S 
m = S 


where 6 


ms 


L^S 


is defined as 



for m ^ S 
for m = S 


This relation follows from the orthogonality properties 



cos ma cos 


ja da = 




+ 2it & 





Pp^(p)pp)( M .)dp = 


l ' — - — 5 

n-j)! 2n+l n,k 


m i J 

m = j / 0 
m = j = 0 


As is explained in the AMLYSIS section, 

2*F 0 q C r ) — q ( r ) j the net neutron current 

the total source strength 2«S Q Q = s Q Q . 


the total neutron flux 
*F^i(r) » and 


Substituting equation (B8) into (B6), multiplying (B6) by 
cos ja P^(p) da dp and integrating over a from -it to at and then 

integrating over p from -1 to 1 (over all directions) and after 
dividing through by 2n 


In this appendix, it is particularly convenient to introduce a 
"modified" delta function notation defined in this way. 
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E £ Mr Cl - * 2)1/2 

n=0 m=*0 ^-1 


jrt 


cos a cos ma cos ja da + 


mF (r) 

/ sin a sin ma cos ja da + 


Na F (r) 
J 


( 5 j,0 + I V>) - S JA (r) ( S J,o + I v) + 


(•* 


Vs,k F j^ Cr) ( 5 j,o + l 6 jjto 


Performing the integration over a, the leakage term in the P3 
approximation can be written 

1 / 1 
1 2j|ti r (l - n 2 ) 1 2 pJ: J 1 (p)P^ m) (lO<¥ E 6 m,p 5 j,l-l 


yy ^ f n + 5 < 


m*P \ 0^P-1 0 


,p+l I m,n r 1 —j 

J P^O 


p=0 


6 5 + 

m,p J ) 1”P 






m F. 


&_ J8, „ , - 5, 


m^n I 


m,pl j,p-l j>P +1 / r 


or 


EEteH^ f ^ fp,» ♦ v-p + 

n=0 b=0 ^-1 ^ J 


n=0 p=0 

oo 3 

EE 

n=0 p=0 


(n-p) 

! 2n+l 

/> 

- n 2 ) p^ p_1 ^(^)P^ p H^) d + i | 

F- + p^\ 

(n+P) 

! 2 

L. P,n I 

(n-p) 

! 2n+l 

/> 

- ii 2 ) 1 pjf +1 V)pjfV) d P < 

[p ' n - P ^ 

(n+p) 

! 2 

1 P> n r 


,p+l 


n=»0 p=0 

where j has positive values only 
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After evaluating the remaining integrals in the leakage term (ref. 
3, pp. 102 to 104) the Boltzman equation can finally be written 


k(k+l) 

4 ( 2k+l ) 


F' - F 

0,k-l 0 


1 6 + 1 5" + F l.k+1 F l,k-l| R 

),k+lj j,l 4(2k+l) |*l,k+l r l,k-l r J 8 j,0 

v i r F , + p F P,k + i pi , _ p F Pi k-i i 

/__j 4(2k+l) p,k+l r p,k-l r 

p=l L J 


+ 


p=0 

N a F 


6 . + 
J,P+1 


2 r 

E p P+l)(k+ P)( F p,k-l - ^ Eli ) ‘ (k-P+l)(k-p)^F^ k+i - 

lk (r) ( 5 j ,0 + 2 6 J? to) = Vs,k F j,k (r) [^0 + I S jjio] + S j, k (r) ( 5 j ,0 + \ 6 J? fo) 


In the approximation, a set of six simultaneous equations in 

terms of F m n (r) are generated by letting j and k independently 
take on the values 0,1, 2, 3 


r 


F m,n( r )= 0 for < 


m or n < 0 
m or n > 3 
m > n 

^either subscript even and the other odd 


Assuming a constant isotropic source of neutrons, the following 
set of simultaneous equations are generated: 


F' . + 

1,1 r 


f l,3 


1,1 


F ' , + - F ' 

1,3 r 1,1 r 


+ 2b J? . = 2S 

0 0,0 “*0,0 


1,1 


^ + 10 b 2 F 0,2 = ° 


2FA r, - 2F1 o + i F ' _ + -2*! + 3b p , = o 


0,0 - " 0,2 


11, 1 


(B9a) 

(B9b) 

(B9c) 
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12FA „ - J FA o - + 7b 3 F i,3 " 0 


0,2 2 2,2 r 


33? ^ 3 

F 3,3 + + 121 


( h ,1 - V) - 2 ( F i,3 - ¥) 


3 .-^i+iob^a 




^ + 7b,F 


3 3,3 


0 


(B9d) 

0 

(B9e) 

(B9f) 


where b^ = Ncr - ®g8g ]^} k = 0, 1, 2, 3. 

By setting F Q 2 and F £ g = 0 in (B9c), equations (B9a) and 

(B9c) reduce to the ordinary cylindrical diffusion equation when they 
are transformed by 


2lt F 0,k _4> 0,k 


F J,k - *j,k 


2jt s o,o “ s o,o 
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APPENDIX C 


CALCULATION OF COEFFICIENTS IN P g APPROXIMATION 
NEUTRON FLUX EQUATIONS 

The term a k appearing in equations (l) is given by 


a k - zA - *s°s,k 


(Cl) 


where the superscripts (which will hereafter be dropped) refer to region 
A. The cr g ^ k appearing in (Cl) is the coefficient of the k th Legendre 

polynomial in the series expansion of the scattering function, that is 


The physical meaning of Og Q and a g ± can be deduced from the 
following considerations: 3 3 

The term N s a s (S2,fl') is the probability per centimeter that a neu- 
tron traveling in the direction S2 ' is scattered into the direction 2 
from following sketch. Scattering is assumed to depend o nl y on the 



By the orthogonality of Legendre polynomials 



(C2) 
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e'- 

er} 

co 

to 


angle 0 (or its cosine) where ft' •£ = cos 0 = p Q . Then a quantity 
N S a sC |i o) d ^O can de d - efined as the probability per centimeter of a 

neutron being scattered from S$' into an elemental ring of area 
2« dp 0 with direction cosines between Pq and Pq + dp^. Then 

_2* 

N s dp Q / or s (^S')cii = %o s (^ 0 )d^ 0 

d|r=o 

or 

N s o s (2,a') = °s(^0^ 

For k = 0, the right side of equation (C2) is the scattering 
function integrated over all velocity directions so that Og q is the 

total microscopic cross section for scattering. 

For k « 1, P-lCm-q) = M- 0 so that Cg^ is a s ^ Q times the average 
cosine p Q of the scattering angle. 

Now 

o s (p 0 ) d ^0 = pC^o^o a s,o 

where a c n is the probability that a neutron is scattered and 

O y U 

p (n Q ) dp Q is the probability that the scattered neutron has the direc- 
tion cosine Pq between Pq and Pq + dpQ . This probability must be 
independent of the coordinate system, therefore 

p(n 0 ) d P 0 - p(v)dv 

where p(v)dv is the probability of scattering into the angle 

2it d(cos 0 C ) in the center of mass system where 0 C is the scattering 

angle in the center of mass system and v = cos 0 C . 

Assumption of spherically symmetrical scattering in the center of 
mass system implies 

r \ I 2ndv 1 

P(v)d-v = — = 2 dV 

Therefore, 

a s(^0^ d|i 0 “ dv 


t 
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Applying the laws of conservation of energy and momentum for elas- 
tic collision of a neutron with a stationary nucleus of mass M, yields 

Vi 

“O * — 175 W 

(ir + 2Mv + 1) 

which relates the cosine of the scattering angle in the laboratory sys- 
tem to that of the center of mass system. Equation (C2) can then be 
written 

1 

P k W dv (C4) 


Equation (C3) enables the integral in (C4) to be evaluated for 
k=0 ,1,2,3 in a straightforward manner with the results 



°S,0 " CT s,o 




V " 3M ° S ,o 


S,2 = °S,0 )q 


5M 

8 


32M 


c“ 2 - u* - (^i) 2 } 


> (C5) 


°S,3 " 0 


J 


The special case of an infinite atomic mass k m (often assumed for 
heavy nuclei) gives Og^ = Og^ 2 = a g g = 0 which follows from equation 

(C5) . Scattering is then isotropic in the laboratory system of reference, 
as well. For light atomic nuclei, the scattering is predominantly for- 
ward in the laboratory system. 


Equations (C5) enable the calculation of a g ^ for the P g approx- 
imation of transport theory. In particular, the values obtained were 
used to calculate a^ (eq. (Cl) ) , b^, and c^. appearing in the flux 

equations for the regions A, B, and C . Only the hydrogen nuclei of the 
water region were assumed to have a finite mass number. 


» 
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TABLE I. - RECTANGULAR GEOMETRY 
[Homogeneous value of R^ = 0.712.] 



Total 

production 

Total 

absorption 

Total fuel 
absorption 

R 

R - 0.712 

Self -shielding factors for 

two -region 

cells 


Transport theory, 






Unbound 

0.1993 

0.1987 

0.1436 

0.7227 

0.0107 

Bound 

.1993 

.1993 

.1449 

.7270 

.0150 

Diffusion theory, P 

u 






Unbound 

.1993 

.2006 

.1438 

.7163 

.0048 

Bound 

.1993 

.1989 

.1432 

.7200 

.0080 

Self -shielding factors for three-region c< 

5 11s 


Transport theory, P^ 






Unbound 

0.1993 

0.2029 

0.1399 

0.6895 

-0.0225 

Bound 

.1993 

.2039 

.1406 

.6896 

-.0224 

Diffusion theory, P^ 






Unbound 

.1993 

.2006 

.1404 

.6999 

-.0121 

Bound 

.1993 

.2026 

.1418 

.6999 

-.0121 
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TABLE II. - CYLINDRICAL GEOMETRY 
[R H = 0.712] 



Production 

per 

(cc)(sec) 

Absorption 

per 

(cc) (sec) 

Fuel 

absorption 

per 

(cc) (sec) 

R 

R - 0.712 

Self -shielding factors for two-region cell 


Transport theory, Pj 






Unbound 

1.000 

0.9880 

0 . 7248 

0.7336 

0.0216 

Bound 

1.000 

.9932 

.7305 

.7355 

.0235 

Diffusion theory, P 2 






Unbound 

1.000 

.9790 

.7047 

.7198 

.0078 

Bound 

1.000 

.9802 

.7070 

.7213 

.0093 
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(b) Three-region. 

Figure 1. - Rectangular cells. 
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(a) Two-region cell. Molecular binding effects neglected. 
Figure 2. - Neutron flux distribution. 
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(b) Two-region cell. Molecular binding effects of water approximated. 
Figure 2. - Continued. Neutron flux distribution. 
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(c) Three-region cell. Molecular binding effects neglected; uranium region shown 
magnified 20 times. 


Figure 2. - Continued. Neutron flux distribution. 
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(d) Three-region cell. Molecular binding effect of water approximated; uranium region shown 
magnified 20 times. 

Figure 2. - Concluded. Neutron flux distribution. 
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(a) Used to obtain equivalent two- 
region cylindrical cell. 



(b) Equivalent two-region cylin- 
drical. cell. 

. • ... * V 

Figure 3. - The cylindrical cell. 
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Figure 4. - Cylindrical coordinate system. 
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Figure 5. - The unbound two-region cylindrical cell having same homogeneous self -shielding factor as slab 
cell (0.712). 
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